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ABSTRACT 

We use a numerical nonlinear multigrid magnetic relaxation technique to investigate 
the generation of current sheets in three-dimensional magnetic flux braiding experiments. 
We are able to catalogue the relaxed nonlinear force-free equilibria resulting from the 
application of deformations to an initially undisturbed region of plasma containing a 
uniform, vertical magnetic field. The deformations are manifested by imposing motions 
on the bounding planes to which the magnetic field is anchored. Once imposed the new 
distribution of magnetic footpoints are then taken to be fixed, so that the rest of the 
plasma must then relax to a new equilibrium configuration. For the class of footpoint 
motions we have examined, we find that singular and nonsingular equilibria can be 
generated. By singular we mean that within the limits imposed by numerical resolution 
we find that there is no convergence to a well-defined equilibrium as the number of grid 
points in the numerical domain is increased. These singular equilibria contain current 
"sheets" of ever-increasing current intensity and decreasing width; they occur when the 
footpoint motions exceed a certain threshold, and must include both twist and shear 
to be effective. On the basis of these results we contend that flux braiding will indeed 
result in significant current generation. We discuss the implications of our results for 
coronal heating. 

Subject headings: MHD, Solar Corona, Current Sheets, Coronal Heating, Flux Braiding 
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1. Introduction 

Accounting for the elevated temperature of the 
corona is a fundamental issue in solar observation 
and theory. For many it is the magnetic field that 
is at the heart of the mystery, acting as the link be- 
tween the convective motions in the photosphere and 
the resultant currents and resistive dissipation in the 
corona itself. Such a link was perhaps first developed 
by Gold (1964), and subsequently formalised to some 
degree by Parker (1972,1983) in his notion of "topo- 
logical dissipation" to describe the process by which 
twisting and braiding of the ambient magnetic field 
can lead to energy release in the corona. Indeed, in 
his 1972 paper, Parker reduced the problem to that 
of considering nothing more than a uniform field con- 
tained between a pair of superconducting plates. This 
latter configuration has become known as "Parker's 
model" , and has proved useful as a testing ground for 
coronal heating theories. 

The currents j produced in the corona must be 
of significant intensity to account for coronal heat- 
ing. This is because the coronal resistivity r\ (typi- 
cally of the order of 10 _12 Mm 2 s _1 ) in the ohmic dis- 
sipation term ijj 2 is so small that only currents of 
a sufficient magnitude can lead to energy deposition 
within the required coronal timescales. To produce 
such currents, the magnetic field must be able to col- 
lapse to sufficiently small scalelengths. Indeed, if we 
take the equilibrium field to contain scalelengths of 
order unity, then we require a collapse by at least 
six orders of magnitude to generate the intense cur- 
rents. Such fine scale current structures we shall refer 
to as "current concentrations" . These are to be con- 
trasted with "current sheets" which are taken to be 
those structures resulting from spatial discontinuities 
in the magnetic field direction, and whose existence 
within the framework of ideal magnetohydrodynamics 
(MHD) is generally considered to be singular, in that 
it is only the presence of finite resistivity that allows 
their spatial extent to be resolved. Current concen- 
trations, on the other hand, can be resolved within 
ideal MHD. However we choose to define these fea- 
tures, it is this vast range of scalelengths they imply 
that certainly makes the task of finding evidence for 
them so challenging. 

From an analytical point of view, Parker's model 
(1972) contends that current sheets will in general be 
created in coronal magnetic fields. This is based on 
Parker's belief that equilibrium is only possible if an 



ignorable coordinate exists in the magnetic field per- 
turbations. Given the random nature of photospheric 
flows, it seems unlikely that the perturbations will 
possess a certain symmetry. Therefore, in order to 
proceed towards an equilibrium, the magnetic field 
must regain this symmetry, and to do so the field 
must alter its topology. Parker suggests that this is 
accomplished through current sheets and magnetic re- 
connection. Since the topology is changed through 
the nonidcalncss of the plasma, the term "topological 
dissipation" was coined. However, van Ballcgooijen 
(1985, 1988a) has pointed out a flaw in Parker's orig- 
inal analysis, and claims in his own right that the 
field can pass through a series of equilibria without 
requiring reconnection. These equilibria arc contin- 
uous and have no symmetry. Hence it is the time 
evolution through these equilibria that can lead to 
coronal heating via current concentrations, and not 
as a result of the coronal field attempting to relax 
to any particular equilibrium state. The jury is still 
out on this issue, however, as evidenced by Parker's 
recent work (1994). 

Numerical attacks have been made on the problem, 
notably van Ballegooijen (1988b, 1990), Mikic et al 
(1989), Longcope and Sudan (1994), Galsgaard and 
Nordlund (1996) and Hendrix and Van Hovcn (1996). 
On the evidence from such modelling it does seem 
plausible that magnetic scalelengths much smaller 
than those proscribed for the imposed flows can result 
in the corona, thereby lending weight to this process 
being a coronal heating mechanism. Furthermore, the 
work of Mikic et al (1989) shows that force-free equi- 
libria without current sheets can be generated for each 
particular distribution of footpoints comprising the 
overall evolution. Longcope and Strauss (1994a, b) 
have examined the formation of current layers with 
a nonzero, but small, thickness. They show, through 
the analysis of of the Jacobian of the ficldlinc map- 
ping, that these current consentrations may easily 
be more than six orders of magnitude smaller than 
the equilibrium length scale. The resulting structures 
they observe bear a marked resembelance to those 
we find in this paper, albeit through the use of dif- 
ferent methods to those presented here. On the ba- 
sis of these few experiments it would seem that the 
weight of evidence is against Parker. However, as we 
have previously noted, the coronal scalelengths actu- 
ally envisaged are so small as to be out of the reach of 
such numerical simulations for the foreseeable future. 
We must therefore tread carefully in extrapolating the 
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numerical results to actual coronal conditions. 

In order to try to alleviate the numerical short- 
fall, resulting from the severe scalelength constraints, 
Craig and Sneyd (1986) (1990) developed a mag- 
netic relaxation technique using a Lagrangian scheme. 
Such a scheme permits the numerical mesh to move 
with the fluid elements and magnetic field lines com- 
prising the plasma. Consequently, as any short scale- 
lengths develop in a region of the plasma, so the 
Lagrangian scheme will aggregate more grid points 
there. This means that for the same number of initial 
numerical nodes as an Eulerian code might have, the 
Lagrangian scheme will be able to resolve much finer 
scales. This is clearly an advantage when the search is 
directed specifically towards finding evidence of cur- 
rent concentrations or sheets. We will therefore be 
using the Lagrangian scheme in this study. This is 
not to say that the Lagrangian code is a panacea for 
all our resolution ills, since we still need to have suf- 
ficient resolution over the whole magnetic structure 
containing any of these current structures. Even al- 
lowing for this latter constraint, there is no doubting 
the efficacy of the Lagrangian scheme in allowing us to 
explore a greater range of parameter space than would 
be accessible to any Eulerian code for the equivalent 
amount of computing power. 

We wish, therefore, to directly test Parker's model 
of current sheet formation using the Lagrangian re- 
laxation scheme. The method allows us to find the 
nonlinear, force-free equilibria consistent with the im- 
posed magnetic footpoint displacements. We are con- 
sequently only interested in the final magnetic con- 
figuration within the constraints of ideal MHD. This 
contrasts with most of the previous numerical stud- 
ies of this problem which were Eulerian and time- 
dependent. It should be noted, however, that al- 
though in the absence of definitive analytic and nu- 
merical diagnostics, we cannot be certain that a so- 
lution with a singular current sheet is necessarily the 
end result of a given numerical experiment, our re- 
sults at least indicate those configurations with the 
potential to develop field singularities. 

The layout of the paper is as follows. In Section 
2 we outline the basic equations and method of so- 
lution, and detail the numerical experiments that we 
have undertaken. Section 3 presents the results and 
analysis from these experiments, and we finish with a 
discussion and some conclusions in Section 4. 



2. Model and Numerical Details 
2.1. Basic Equations 

For this particular study we adopt the low-/? ap- 
proximation for the solar corona and look for equilib- 
rium solutions of the equation, 

(VxB)xB = 0, (1) 

where B is the magnetic field. These equations are 
nonlinear. To find solutions of these equations we 
follow the relaxation method detailed by Craig and 
Sneyd (1986) (1990) and use a momentum equation 
which ignores inertial effects in favour of a dominant 
frictional term proportional to the fluid velocity v, 
and takes gas pressure to be negligible in line with 
the low-/? approximation. This simply leaves us with, 

v=jxB, (2) 

where the current j = V x B. In this form the momen- 
tum equation guarantees relaxation toward a state of 
lower magnetic energy. The equation set is completed 
by noting that B/p satisfies, 

^2 = (5.v,v, ,3, 

Dt p p 

where p is the fluid density, D/Dt is the so-called 
convective derivative, and 

V • B = 0. (4) 

The details for the solution of these equations us- 
ing Lagrangian coordinates can be found in Craig 
and Sneyd (1986) (1990). The vital step is the re- 
alisation that the whole problem can be formulated 
in terms of the Lagrangian position of a fluid ele- 
ment x, say, by replacing v by Dx/Dt, and using 
the so-called Cauchy solution for equation (|j) which 
expresses B/p in terms of the equilibrium magnetic 
field Bo, density po, and x itself. The resultant sys- 
tem of equations is parabolic, as shown by Craig and 
Sneyd (1990). As a consequence, unconditionally sta- 
ble numerical schemes must be used to solve the sys- 
tem if inefficiently small time steps are to be avoided. 
The implicit formulation devised by Craig and Sneyd 
(1986) gives the requisite numerical stability, and can 
be solved using alternating direction implicit (ADI) 
techniques. It should be noted that although La- 
grangian, the equations themselves are solved on a 
fixed grid using centered differences. 
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As for all classical methods, the number of itera- 
tions required for convergence scales as some power 
of the number of points (N) in the system, typically 
as N 2 . Thus for numerical schemes in three dimen- 
sions, doubling the size of the computational domain 
requires approximately 64 times the number of iter- 
ations. It was found that for this simulation grids of 
over 32 3 required an excessively long computational 
time. In order to somewhat alleviate this problem we 
have implemented a nonlinear multigrid algorithm for 
grid sizes of 32 3 or more. For such a scheme the num- 
ber of iterations scale at a fraction of that of classical 
methods thus allowing high resolution calculations to 
be conducted in realistic times. Our method is out- 
lined in the appendix. 

2.2. The Numerical Experiments 

Our equilibrium configuration is exactly that en- 
visaged by Parker (1972), i.e., a uniform, vertical 
magnetic field Bo = B^z, with its ends anchored 
into superconducting plates at z = ±L Z . The ini- 
tial positions of the fluid elements are distributed uni- 
formly throughout a domain bounded by the super- 
conducting plates in z, and by distances x = ±L X and 
y = ±L y . In the previous experiment of the nonlinear 
evolution of the kink instability detailed by Craig and 
Sneyd (1990), the fluid elements lying in the bound- 
ing planes in x and y were taken to be fixed. There- 
fore only fluid elements inside the computational do- 
main were free to move. In our initial experiments 
we adopted the same boundary conditions. However, 
the nature of our experiments meant that we gener- 
ated boundary layers at the x and y bounding planes 
which tended to mask the physics of importance. To 
overcome this problem, the x and y directions are 
taken to be periodic in the displacements of the fluid 
elements, leaving the only fluid elements fixed to be 
those in the planes z — ±L Z . Apart from the ini- 
tial equilibrium, this is the only major change to the 
classical Lagrangian code from its previous incarna- 
tion under Craig and Sneyd (1990). However for the 
systems with large shears and large numbers of grid 
points a nonlinear multigrid scheme has also been 
implemented (sec Longbottom, Fielder and Rickard 
(1997) for further details). 

Starting from the equilibrium distribution just de- 
tailed, we impose displacements of the fluid elements 
within the planes z — ±L Z . These displacements take 
the form of shears of equal magnitude but opposite 
direction on the two driving boundaries. Once dis- 



placed, these elements are then held fixed, and the 
remaining fluid elements are allowed to relax towards 
a new equilibrium configuration consistent with these 
new boundary conditions. We say that we have a 
converged solution when both our norm of the resid- 
ual force has fallen to at least O(10~ 3 ) (starting at 
O(10)) and the maximum current has converged to 
two decimal places. By solving the ideal MHD equa- 
tions, we guarantee that the magnetic field lines are 
"frozen-in" to the fluid. The equilibrium field lines in 
this case are vertical lines running from the bottom 
to the top boundary in z, so that fluid elements ini- 
tially sharing the same values of x and y lie on the 
same field line. As we have intimated, they will lie on 
the same field line throughout the relaxation. We can 
therefore plot any relaxed field line simply by drawing 
a line through that same set of fluid elements. This 
results from our choice of equilibrium, and using the 
Lagrangian coordinate system. 

For every single solution, we attempt to detail nu- 
merical convergence in the solution by repeating the 
same relaxation with ever increasing numbers of fluid 
elements. In this way we hope to separate out equi- 
libria that are non-singular (smooth) from those that 
are singular, i.e., containing current sheets. If the 
equilibrium is singular, then increasing the numbers 
of fluid elements will not result in a convergent max- 
imum current in the current sheet. The scalelength 
of such a singular feature will always be smaller than 
that attainable numerically, and so the current locally 
will continue to increase. However, the total current 
in the computational domain will converge, showing 
that the global solution is convergent, apart from the 
locality of the singularity. On the other hand, the 
nonsingular solutions will converge to a smoothly re- 
solved structure throughout the whole domain. Even 
with the Lagrangian code, the ultimate lack of nu- 
merical resolution means that the convergence prop- 
erties of any equilibrium solution are all that we can 
base our physical results upon. The relaxation ex- 
periments described by Billinghurst et al (1993) to 
examine current sheet formation in two-dimensional 
fields provide a substantive example of this process, 
and it is indeed their convergence tests that we have 
employed here. It should be noted that although in- 
creasing the number of grid points reduces the nu- 
merical error for a particular solution, the solution 
itself depends on the position of the Lagrangian grid 
points. Within the confines of one dimensional the- 
ory it can be shown that the motion of the Lagrangian 
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grid allows a better representaion of the true solution 
as compared to an Eulerian grid with the equivalent 
number of grid points. The maximum current that 
may be resolved scales as the square of the number 
of grid points. The behaviour in a higher number of 
dimensions is not so obvious, but for many cases an 
improved representation of the true solution over the 
equivalent Eulerian code is found. 

Having outlined the basis for our numerical exper- 
iments, we now need to consider possible footpoint 
motions. The key ingredient is to produce scalelength 
collapse of the magnetic field within the volume as 
a result of smooth displacements on the boundary. 
This is the crux of Parker's model (1972). Clearly, 
fine scale motion on the boundary will manifest itself 
as fine scales within the volume. However, the ob- 
servations suggest that the scale of convective motion 
in the photosphere is much greater than we envisage 
for the current sheets. It must therefore be left up to 
the coronal magnetic field to produce the short scale- 
lengths of its own accord. Parker (1972) believed that 
the coronal magnetic field always does, and hence the 
notion of topological dissipation. 

Our particular choice of footpoint motions are 
motivated by the work of Galsgaard and Nordlund 
(1996). They conclude that the most important fac- 
tor leading to current sheet formation is one shear in 
one direction followed by another in the orthogonal di- 
rection. This process is deemed to be sufficient to pro- 
duce the exponential current growth predicted by van 
Ballegooijen (1986) on the basis of a random series of 
shearing motions. As Galsgaard and Nordlund (1996) 
note, it is the underlying field and flow topology re- 
sulting from the shear plus a shear that inevitably 
leads to the exponential growth. It seems clear that 
the alternating shear flow profiles advocated by van 
Ballegooijen (1988a), and tested by van Ballegooijen 
(1988b) and Mikic et al (1989), will still lead to the 
exponential growth, since the sequence of footpoint 
motions is indeed composed entirely of the basic ele- 
ment, i.e., a shear plus a shear. To simplify the prob- 
lem as much as possible, then, we choose to examine 
the possible effect on the coronal field of a shear in 
one direction, followed by a shear in the same, or in 
another, direction. It is these basic elements that we 
are going to focus on here. 

In the numerical experiments, distances are nor- 
malised to the initial length of each field line, and 
we take the normalising magnetic field strength and 
plasma density to be those in the equilibrium. For 



the experiments we have conducted to date, we use 
the unit cube as the initial volume of our numerical 
domain, and we take the field strength to be one and 
density to be 0.05. These provide the reference points 
when we later detail results from the relaxed equilib- 
ria. 

To maximise the spatial resolution available to us, 
we will simply examine the equilibria associated with 
the footpoint displacements shown in Figures Al(a)- 
(f). Figure Al(a) shows the initial equilibrium foot- 
points. At every grid point (at the intersection of grid 
lines) lies the footpoint of each magnetic field line. 
The Lagrangian code follows these particular points, 
so that after each successive distortion of the grid we 
can exactly identify whe re t he original footpoint lo- 
cations shown in Figure Al (a) have been moved to. 
Figure |Al| (b) results from a shear of magnitude 0.8 
appli ed parallel to the y-axis. The remaining Fig- 
ures AL(c)-(f) result from shears of magnitudes 0.1, 



0.3, 0.5, and 0.7, respectively, applied parallel to the 
x-axis after having previously applied the shear re- 
sulting in Figure ^l|(b). We are therefore examin- 
ing the effect of shearing in one direction, immedi- 
ately followed by shearing in the perpendicular di- 
rection. We hope to show that this is sufficient to 
produce scalelength collapse within the plasma vol- 
ume. It should be noted that the displacements on 
the boundaries L = ±L Z are of the same magnitude 
but of opposite sign. 

From the evidence of some of the gross grid distor- 
tions shown in the previous Figures, the reader may 
be surprised that we can claim to obtain relaxed equi- 
libria from such complex arrangements. We should 
point out, however, that we have shown the physi- 
cal grid displacements, and not the uniform grid on 
which the resulting Lagrangian equations are actually 
solved. As we shall show, the more distorted the sys- 
tem, the more difficult the relaxed equilibrium is to 
obtain. Nevertheless, given sufficiently robust relax- 
ation techniques, the Lagrangian system should allow 
us to approach the equilibrium associated with each 
set of footpoints arbitrarily closely. 

Although we have detailed the footpoints, there 
is still the question about what to use as the start- 
ing points for the fluid elements within the compu- 
tational volume. As Craig and Sneyd (1990) have 
intimated, we have a nonlinear equation to solve, and 
it is therefore possible that differing initial conditions 
for the same footpoint distribution will result in dif- 
ferent relaxed states. While it would be interesting to 
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explore this hypothesis, we are faced with the more 
practical problem of obtaining any relaxed equilib- 
ria at all. This is because not all initial conditions 
lead to solutions, despite the unconditional numerical 
stability of the Lagrangian code. It seems that some- 
times the code finds an initial condition too complex 
to be able to sensibly unscramble, and fails. The best 
strategy we have found to date is to use as an initial 
guess the relaxed solution from the previous state in 
the sequence of footpoint displacements, with the new 
footpoint displacements prescribed on the boundaries 
z = ±L Z . So, although we view each equilibrium as 
an entirely separate solution to all the others, there 
is a sense in which they can be viewed as a sequence 
connected together. The "dynamical" experiment of 
Mikic et al (1989) falls into this category. 

Having applied the footpoint displacements to the 
end planes at z = ±L Z , these are then held fixed, 
and the remaining fluid elements are able to move in 
such a way as to approach a new equilibrium consis- 
tent with the footpoint distribution. How close we 
are to the new equilibrium is measured by the maxi- 
mum Lorentz force j x B within the domain. Having 
obtained a "suitably good" solution, we then calcu- 
late the maximum current within the domain, and the 
total integrated current strength |j| over the whole 
domain. We then repeat the whole experiment for in- 
creasing numbers of fluid elements within the domain. 
The whole execise is then repeated for each new dis- 
tribution of footpoints. Based on these results, we are 
then in a position to comment on the smoothness or 
otherwise of a given equilibrium solution. 

Our ability to catalogue each equilibrium accu- 
rately is limited firstly by the available computing 
power, and secondly (and perhaps more importantly) 
by the speed at which the equilibrium solution is 
approached numerically. The last feature is com- 
pounded by increasing numbers of fluid elements and 
the increasing complexity of each equilibrium solu- 
tion. This means that, even with the implementation 
of multigrid techniques, the equilibria associated with 
the most distorted grids are not only the hardest to 
solve computationally, but also, because of their in- 
herent problems, are the ones on which we have the 
least amount of numerical convergence. 

3. Results 

The summary of all our results is encapsulated in 



of the current J m obtained within the computational 
domain (always at the centre of the mid-plane in z) 
in the best relaxed solution, against a measure of the 
number of grid points in the simulation n a . Here the 
total number of points in the simulation is equal to 
n^, and the measure of the best relaxed solution is 
that with the lowest maximum value of the Lorentz 
force j x B obtained so far. The number labelling each 
curve shows the amplitude of the shear applied in the 
^-directio n, w hile the unlabcllcd curve refers to that 
in Figure Al(b). Converged solutions (either using 



Figure A2 . This plots the maximum absolute value 



the classical relaxation method or nonlinear multi- 
grid) are plotted with an *, those that are under es- 
timates of the current are plotted with a A. 

Figure A2 reveals that j m is tending to saturate for 
the lower set of shears applied, i.e., up to and includ- 
ing 0.5. We would therefore conclude that these reveal 
smooth equilibria. However, for the shears labelled by 
0.6 and 0.7 the value of J m is tending to increase at 
least linearly with n a , so that within the confines of 
our resolution we would conclude that these point to 
the possibility of (at least) significant current concen- 
trations, if not current sheets themselves. It should be 
noted that, within the confines of a one dimensional 
Lagrangian representation, a rigorous condition for a 
truly singular current sheet may be obtained. This 
states that the maximum current J m must scale as 
n^. If this is taken as a generally necessary condition 
for current sheet formation, then only equilibria with 
a shear of 0.7 exhibit singular behaviour, J m scaling 
approximately linearly with n a for a shear of 0.6. It 
is not however obvious how the effect of working in a 
higher number of dimensions would modify the con- 
dition. Thus we refer to both shears of 0.6 and 0.7 as 
divergent solutions. 

Although we have a broad range of footpoint dis- 
tributions, the generic forms of the relaxed solutions 
we see divide themselves into two types, at least with 
respect to the structure of the current. The first is as- 
sociated with shearing only in the y-direction. Here 
there is no specific peaking of the current distribution. 
Rather we have a relatively broad curr ent structure 
filling most of the domain. As Figure A2 shows, we re- 
quire only relatively modest numbers of fluid elements 
before we confidently predict that the equilibrium is 
a convergent one with respect to n a . Indeed, for cer- 
tain parameters, this equilibrium may be calculated 
analytically (Hood, priv. comm.). For these cases 
the numerical method reproduces the true solution to 
within the given tolerance. The basic field structure 
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as a function of x of these solutions is shown in Fig- 
ure |Al| Figure [A^(a) shows the variation of B x with 
the height z. The bold curve with the largest ampli- 
tude is B x close to the lower z boundary. As you go 
up in z the amplitude of B x gradually falls, resulting 
in the second bold curve plotted at about half way 
between the lower z boundary and the mid-plane. In 
the mid-plane, B x is zero. From the mid-plane up- 
wards the variation in B x is a mirror image of that 
below the mid-plane, with the sign of B x reversed (as 
indicated by the dashed curves). There is only a sin- 
gle curve shown for B y in Figure [A^(b) as the relaxed 
By is practically uniform in z. The variation in B z is 
similar to that for B x , i.e., a weak boundary layer at 
the top and the bottom, with hardly any z-variation 
about the mid-plane. The only difference is that there 
is no sign change in B z , hence the two curves, the 
larger amplitude one in the mid-plane, and the other 
close to the z-boundaries. 

We also note that the variation in x of B y exactly 
matches of that of the imposed footpoint displace- 
ment, i.e., proportional to sin(27ra;), whereas B z has 
half the wavelength. This variation in B z results from 
the need to balance the excess magnetic pressure in 
the y-direction produced by the footpoint shear. B z 
therefore has to match the x- gradient of B y 2 , result- 
ing in the x-profile we see. To produce this B z profile 
the fluid elements have to displace themselves in the 
x-direction in order to generate regions of either in- 
creased or decreased B z . Since our fluid elements are 
tied to the field lines, this x-motion naturally results 
in a B x that has the same wavelength as B z , but is 
7r/2 out of phase. 

These first type of solutions are therefore almost 
wholly dependent on x alone about the mid-plane, 
with most of the variation in z occurring towards the 
top and bottom planes. The lowest energy state for 
our boundary conditions seems to be one in which an 
almost one-dimensional structure extends over most 
of the domain, with most of the stress taken up in the 
weak boundary layers. This is reminiscent of results 
reported by Browning and Hood (1989) and references 
therein in the context of twisted flux tubes in which, 
as here, the equilibrium is practically one-dimensional 
everywhere, except for boundary layers near the line- 
tied boundaries. This seems to be the lowest energy 
configuration with respect to one in which, as might 
have been anticipated, there is a gradual variation 
in the properties of the equilibrium over the whole 
length of the structure. 



The second (and most important) type of current 
structure is revealed in all the other relaxed solutions, 
once a second shear in the x direction has been ap- 
plied. The only thing that distinguishes each of these 
particular solutions is the amplitude of the current ob- 
tained. Each solution contains a strong current con- 
centration restricted to a narrow slice passing through 
the centre of each z-plane, which then rotates slowly 
as increasingly higher z-planes are passed through. 
The amplitude of the maximum current in each plane 
in this feature is practically constant along its length, 
peaking in the central plane at z = 0, and diminish- 
ing only very gradually towards the two ends. The 
main component of the current is in the z-dircction. 



A representation of this feature is shown in Figure A4 



which shows the isosurface at 50 per cent of the max- 
imum current in j z . Larger amplitudes of j z are en- 
closed within this isosurface, and this surface there- 
fore reveals the extent of the current concentration 
(and consequently its halfwidth, that we define later). 
The effect of the shear is apparent in the twist of the 
structure. We also note that there is just a single 
isosurface at this level in the volume. 

Representative isosurfaces in the same relaxed so- 
lution for the currents j x and j y are shown in Fig- 
ure A5(a) and Figure |A5| (b) , respectively. Comparing 
with Figure A4 we see that the isosurfaces for j x and 
j y exhibit the same twist, but that instead of a single 
isosurface we now see the presence of two isosurfaces 
that, to a first approximation, tend to wholly enclose 
the 50 per cent isosurface for j z . In common with j z , 
however, j y shows a fairly uniform spatial distribution 
in z, whereas j x has a broader distribution near the 
bounding planes in z, which then narrows and rotates 
as the mid-plane is approached from either side. 

The explanation for the second type of current 
structure is fairly simple, and has been detailed by 
Galsgaard and Nordlund (1996). It all boils down to 
how the field lines wrap themselves around one an- 
other as a result of the two sets of shears, and the 
fact that they are line-tied at each end. A schematic 
of t he r esult of wrapping the field lines is given in Fig- 
ure A6, which shows the two planes to which the foot- 



point displacements have been applied, along with the 
resultant distribution of field lines (labelled with ar- 
rows). The field line about which all the others wrap 
is shown as the vertical dashed line. The bold field 
lines lie wholly in front of the vertical line, while the 
dotted field lines lie wholly behind both these and the 
vertical field line. The result of the magnetic tension 
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and the line- tying forces the field lines to compress on 
top of one another, producing the strong currents that 
we see. In particular we note that despite the shear 
the main component remains that in the z-direction. 
Hence the single isosurface encompassing the current 
layer shown in Figure A4. However, the components 
j x and j v retain information relating to the imposed 
shears, and hence have a change of sign on passing 
through the current layer in j z . Therefore, j x and 
j y have local maxima on both sides of the j z current 
layer, and this explains the fact that they each reveal 
isosurfaces either side of this layer. 

While the isosurfaces reveal the general structure 
of the equilibrium current, we have lost a sense of 
direction in them. This can be recovered by consider- 
ing some specific contours. In particular, Figure A7 
shows contours of j x in planes of constant z starting 
from the bottom plane, and working up to the mid- 
plane. There is no need to go further in z since, as 
the isosurfaces show, the structures simply continue 
to rotate from those below the mid-plane, albeit in the 
opposite direction. The dotted contours are for neg- 
ative current values. We now clearly see the effect of 
shear, with j x reversing direction about a specific line 
below the mid-plane, and actually reversing about the 
line y = in the mid-plane itself. For j x , then, the 
structure in z is tantamount to taking the mid-plane 
contours and slowly "sliding" them past one another 
as you descend or ascend through the planes. 

The equivalent contours for j y are shown in Fig- 
ure A8. The mid-plane contours now reverse across 
the line x — 0. As you move through the planes in 
z the contours rotate together, unlike the "sliding" 
observed in j x . The sets of contours for j x and j y 
not only have a certain spatial symmetry, but also 
their amplitudes in that in each contour the maxi- 
mum (positive) contour is always the same magni- 
tude as the minimum (negative) contour. The abso- 
lute maxima are (a) 16.4, (b) 9.9, and (c) 5.3 for j x 
in Figure O and (a) 7.8, (b) 10.2, and (c) 12.0 for j y 
in Figure AS . We see that j x is tending to dominate 
j y close to the lower boundary, while the opposite is 
true by the time we reach the mid-plane. 

As we would anticipate from the isosurface plot, 
the contours of j z shown in Figures |X^(a)-(c) are 
dominated by the vertical (positive) component. We 
also recover the sense of a single feature that rigidly 
rotates as we progress from the lower to the upper 
boundary. Specifically, the maximum and minimum 
contour values are (a) 36.5 and -9.9, (b) 37.8 and - 



4.6, and (c) 37.5 and -3.5, respectively. These con- 
tours also reveal how remarkably uniform the current 
concentration in j z is with regard to its length and 
breadth in each z-plane. We can therefore simply use 
the dimensions of the current concentration in the 
mid-plane as a reference point. Furthermore, since we 
obtain equilibrium solutions, then the current j must 
be everywhere parallel to the magnetic field B. Hence 
knowing the field line topology immediately gives us 
the current topology. 

As we have indicated, the main aim of these exper- 
iments is to determine whether or not current concen- 
trations or sheets are possible in these sheared equilib- 
ria. Since the anticipated scalelengths are well below 
those accessible numerically, it is only through tests 
of the numerical convergence of our solutions that we 
can make predictions about the behaviour of the ac- 
tual coronal fields. It is therefore important to have 
suitable reference points in each solution that can be 
reliably measured, so as to allow the checks for con- 
vergence. Global checks of the total current are one 
measure, as well as any local current maxima. Here, 
the mid-plane solution provides such a useful refer- 
ence point. In particular, the current concentration 
associated with j z is clearly produced by the shorten- 
ing of the scales in the x direction. We will therefore 
use the distance across the line x = between the two 
points where the current falls to half its peak value 
along the line y = as a relevant indicator of conver- 
gence. 

Having defined this half-width length scale, we can 
map this on to the the number of grid points that an 
equivalent Eulerian code would require to gain the 
same resolution as the Lagrangian code. These are 
shown in Table |l|. Obviously one could argue that 
an Eulerian code could used a number of methods, 
for example a stretched grid, to gain more resolution. 
However an equivalent method could normally be ap- 
plied to the Lagrangian code. For the largest shears 
we have been able, with the use of the Lagrangian for- 
mulation coupled with nonlinear multigrid, to resolve 
current structures of approximatley l/1000t/i the size 
of the computational domain. This is an order of 
magnitude better than an equivalent Eulerian code 
(see for example Galsgaard and Nordlund (1996)). It 
should be noted, however, that the Lagrangian for- 
mulation by definition can only examine ideal MHD, 
and unlike resistive Eulerian codes nothing can be 
said about reconnection or subsequent evolution once 
resistive effects come into play. The resulting high 
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resolution can however point to whether the current 
structures seen in resistive codes are indeed "current 
sheets" or just current concentrations that would al- 
low slow reconnection at solar values of resistivity. 
Here slow is with regard to observed timescales. 



Figure AlC shows a plot of maximum current 
against half-width (A) as the number of grid points 
is increased. Only the two shears that appear to give 
a divergent current are shown. A tentative extrap- 
olation may be carried out to the width of a struc- 
ture allowed by solar resistivity. Fitting a power 



law to the data in figure A10 gives the relationship 
3 m — 3.78/A 81 . This corresponds, for a lengthscale 
collapse of about six orders of magnitude, to a cur- 
rent of the order of 3 x 10 5 B Q /(^L ) Am~ 2 , where B 
and Lq are characteristic magnetic field strengths and 
field line lengths. Taking B = 10G and L„ = 10 6 m 
this gives a current density of 2 x 10 6 Am -2 . If the 
energy within the current sheet is released by Joule 
dissipation then this total energy may be calculated 
as 



3D = 



j 2 1 a dv 



(5) 



giving jo = 4 x 10 18 J. For dissipation by reconnec- 
tion acting on the Alfven time scale (say the order 
of 100 seconds) we would get an energy flux of ap- 
proximately 4 x 10 4 Js _1 m -2 over the area of the 
box. Given the number of assumptions in making the 
above calculation, it is still interesting to see that this 
value is comparable with energy flux observed in ac- 
tive regions of 10 Js m (see for example Withbroe 
and Noyes (1977)). 

A further measure of how close our solutions are 
to the true equilibrium can be gained by plotting the 
integral of | j | over whole the box and the magnetic en- 
ergy within the box against the resolution used. Both 
of these quantities should converge with increasing 
resolution for an equilibrium solution, whether a cur- 
rent sheet is formed or not. We show these plots in 



Figure All. It can be seen, in all cases, that there is 
a smooth convergence with increasing number of grid 
points. 

4. Discussion and Conclusions 

In this paper we have addressed the issue of the 
production of current concentrations or sheets in non- 
linear, force-free equilibria in coronal magnetic fields. 
Previous work suggests that the combination of a 



shear plus a shear is sufficient to do so, and we there- 
fore focus on such footpoint displacements. Our re- 
sults to date show that for certain amplitudes of the 
second shear, convergent, well-resolved equilibria can 
be generated. However, beyond a critical point it 
seems that the equilibria are divergent, in the sense 
that for increasing numbers of fluid elements we find 
convergence in the total current and magnetic en- 
ergy in the system, but that a local current peak 
(and its associated scalelength) are not convergent. 
On extrapolating these latter results to scalelengths 
expected in the coronal environment, we would pre- 
dict that currents of sufficient magnitude to account 
for significant heating can be generated. Whether 
these equilibria are singular in the sense of leading 
to current sheets we cannot determine. Nevertheless, 
the current anticipated is sufficiently concentrated to 
make the point academic; layers of resistively heated 
plasma will be present. 

Furthermore, by use of a Lagrangian approach, 
we have been able to access scalelengths well below 
those attained by any Eulerian code to date. The La- 
grangian code allows fluid elements to migrate into 
those regions where scalelengths are reducing, and 
therefore makes it particularly suitable for such equi- 
librium experiments. Field lines are also simply plot- 
ted by tracing through the positions of the fluid ele- 
ments; there is no recourse to the use of integration 
routines etc. One disadvantage of the classical relax- 
ation method is its slow rate of convergence for the 
most intensive solutions (in terms of complexity and 
number of grid points). Through the use of nonlin- 
ear multigrid methods we have been able to partially 
overcome this problem. Thus converged solutions to 
these most important experiments can be calculated 
in finite times. 

If the reader is seeking a black and white answer 
to whether Parker or van Ballegooijen is correct on 
the question of the nature of coronal equilibria, then 
this is not the place to look. There is certainly no 
degree of symmetry assumed in our simulations; they 
are fully three-dimensional. Nevertheless our imposed 
footpoint distributions have a certain symmetry to 
them. As for the equilibria, we would contend that 
we find both well resolved, fully three-dimensional so- 
lutions, as well as those that demonstrate a level of di- 
vergence. This puts us firmly on the fence. Whatever 
one's views, we do present compelling evidence for the 
generation of significant current concentrations in the 
corona. 
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The question of the timescale required to produce 
such structures is an important one. We only see the 
final product in terms of the relaxed solution once 
the footpoint displacements have been imposed. It 
must be left to others to determine whether a fully 
dynamical corona can do likewise. Nevertheless, as 
pointed out by Galsgaard and Nordlund (1996), the 
timescale for the footpoint motions is a least a factor 
of a hundred less than that of the characteristic Alfven 
transit time, and one would therefore expect a signif- 
icant amount of time for near equilibrium conditions 
to evolve. Future dynamical experiments with such 
separation of timescales will have to be performed to 
demonstrate the truth or otherwise of this assertion. 

As ever, we have only begun to scratch the sur- 
face. Our bias towards the Lagrangian approach en- 
sures that we will continue to employ such methods 
time dependently, hopefully to answer some of the 
questions of timescale. We are also following up the 
work of van Ballegooijen (1988b) with an Eulerian ap- 
proach using the Euler potentials. With the advent of 
nonlinear multigrid methods and more powerful work- 
stations, it seems entirely appropriate to see how far 
we can push the Eulerian codes in this context! 

The authors wish to express their gratitude to 
Klaus Galsgaard for his prompting on the issue of 
flux-braiding and to the referee for a number of use- 
ful comments. 



A. The Implementation of Multigrid 

Here we outline our implementation of nonlinear 
multigrid with respect to the above problem. More 
detailed accounts of multigrid can be found in Brandt 
(1977), McCormick (1987) and for M.H.S. relaxation 
theory Cally (1991), Fiedler (1992), Longbottom, 
Fiedler and Rickard (1997). 

The key to understanding multigrid methods is 
to know why classical methods fail. Classical relax- 
ation methods are extremely efficient at smoothing 
out short wavelength errors, those on the scale of the 
grid spacing. However they are poor at smoothing 
longer wavelength errors. The trick with multigrid 
is to restrict the problem on to successively coarser 
grids, smoothing the errors in each case, and then 
prolongate the more accurate solution back on to the 
finer grids. As smoothing on the coarser grids takes 
significantly less time than smoothing on the finest 



grid such a multigrid sweep, smoothing all wavelength 
errors, is much more efficient than a number of iter- 
ations on the finest grid smoothing only short wave- 
length errors. 

We implement a full approximation scheme (F.A.S.) 
fixed schedule multigrid, including F-cycles. This pro- 
ceeds as follows: 

For solving F(x) = j x B = 0, we start with a first 
guess for x (x F ) and (pre) smooth on the fine grid 
to find the next best estimate x F +1 . We then com- 
pute the residual rp — F(x F +1 ), and restrict both 
x^, +1 and rp to the next coarsest grid by a second 
order three dimensional weighting Yc — 1Z(yf) and 
x.q = lZ(x F +1 ). Here 1Z represents the restriction 
operator. The approximate solution to the equation 
F(x c ) = r c + F(x%), (x£ +1 ) is then calculated on 
this coarser grid (by recursive multigrid calls to suc- 
cessively coarser grids). We then prolongate the cor- 
rection to the solution, Xp +1 — Xp, back up to the finer 
grid by trilinear interpolation and add it to the best 
estimate on the fine grid x F +2 = x" +1 +7>(x£ +1 -xg,). 
Finally a further (post) smooth is carried out using 
x^, +2 as the initial guess for the relaxation. This com- 
pletes one multigrid sweep. 

In general the efficience of the process is dependent 
on the order in which the coarser grids are visited. 
We find, for this problem, that an F-cycle is the best 
choice. Thus for a four leveled system (labeled 1 to 
4 from coarse to fine) one multigrid F-cycle visits the 
grids in the order 4-3-2-1-2-1-2-3-2-1-2-3-4. 

Ideally for Multigrid the number of iterations should 
scale independently from N (the number of points on 
the finest grid), however in many complex 3D non- 
linear cases the scaling is with some power of N that 
lies between and 2. We find here that the number 
of iterations scale like N as opposed to N 2 for the 
classical method. 



Figure A12 compares the convergence of the clas- 
sical and multigrid versions of the ADI relaxation 
method. The residual force measures the error in the 
solution over the whole domain. One work unit (WU 
is the cpu time for one iteration on the finest grid. 
The examples shown are for grids of 33 3 and 65 3 
points after the second (x-axis) shear of magnitude 
0.2 has been applied. In both cases the convergence 
of the classical method soon saturates (once short 
wavelength errors have been smoothed), however the 
multigrid method continues to converge. Also the 
scalings with number of grid points for both classi- 
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cal and multigrid methods can be seen. 



Lagrangian Effective Eularian resolution for given shear: 



resolution 


0.8 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


21 3 


26 3 


39 3 


42 3 


46 3 


48 3 


50 3 


52 3 


31 3 


26 3 


46 3 


58 3 


70 3 


82 3 


92 3 


96 3 


33 3 * 


26 3 


47 3 


60 3 


76 3 


85 3 


118 3 


120 3 




26 3 


50 3 


68 3 


92 3 


134 3 


200 3 


308 3 


51 3 


27 3 


54 3 


77 3 


96 3 


188 3 


286 3 


556 3 


61 3 


27 3 


58 3 


85 3 


133 3 


234 3 


462 3 




65 3 * 


27 3 


58 3 


86 3 


136 3 


240 3 


521 3 


943 3 



Table 1 : The effective number of grid points needed by an Eulerian code to resolve the current structures obtained 
using the Lagrangian code. Resolutions marked * were calculated using nonlinear multigrid. 
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Fig. Al. — Footpoint displacements applied on the 
plane z = —L z . The displacements on the plane z = 
+L Z are the mirror images of those shown, reflected 
in the planes x = and y = 0. 

Fig. A2. — Maximum current j m obtained in the best 
relaxed solution for each distribution of footpoints, 
plotted against n a (see text). 

Fig. A3. — Magnetic field structure for the first type 
of relaxed solution. The magnetic fields shown are (a) 
B x , (b) By, and (c) B z . 

Fig. A4. — Isosurface of the j z current component at 
50 per cent of the maximum current in the equilibrium 
solution. The vertical direction is z going between 
—0.5 and 0.5. The remaining directions are x and y 
forming a right-handed set and varying between —0.3 
and 0.3. 

Fig. A5. — Isosurfaces of the (a) j x , and (b) j y cur- 
rents at 25 and 42 per cent of the maximum currents 
for each component, respectively. The cube orienta- 
tion and dimensions are as the previous Figure. 

Fig. A6. — Schematic of field-line wrapping resulting 
from the application of two shears. The field lines are 
labelled with an arrow. 

Fig. A7. — Contours of the current j x in the planes 
(a) z = -0.5, (b) z = -0.25, and (c) z = 0. 

Fig. A8. — Contours of the current j y in the planes 
(a) z = -0.5, (b) z = -0.25, and (c) z = 0. 

Fig. A9. — Contours of the current j z in the planes 
(a) z = -0.5, (b) z = -0.25, and (c) z = 0. 

Fig. A10. — Maximum current against half-width 
for a varying number of grid points. Only the two 
second shears, 0.6 (A) and 0.7 (x) that have diverging 
currents are shown. 

Fig. All. — The convergence of the total integrated 
current and magnetic energy as the number of grid 
points are increased. 

Fig. A12. — Residual force against WU for a shear 
of 0.2 with 32 3 (*) and 64 3 (x) intervals, respectively. 
Dashed curves show the classical method, solid curves 
nonlinear multigrid. 
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